Hybrid motion artifact detection and correction approach for functional near-infrared spectroscopy measurements

Abstract. Significance: Functional near-infrared spectroscopy (fNIRS) is a promising optical neuroimaging technique, measuring the hemodynamic signals from the cortex. However, improving signal quality and reducing artifacts arising from oscillation and baseline shift (BS) are still challenging up to now for fNIRS applications. Aim: Considering the advantages and weaknesses of the different algorithms to reduce the artifact effect in fNIRS signals, we propose a hybrid artifact detection and correction approach. Approach: First, distinct artifact detection was realized through an fNIRS detection strategy. Then the artifacts were divided into three categories: BS, slight oscillation, and severe oscillation. A comprehensive correction was applied through three main steps: severe artifact correction by cubic spline interpolation, BS removal by spline interpolation, and slight oscillation reduction by dual-threshold wavelet-based method. Results: Using fNIRS data acquired during whole night sleep monitoring, we compared the performance of our approach with existing algorithms in signal-to-noise ratio (SNR) and Pearson’s correlation coefficient (R). We found that the proposed method showed improvements in performance in SNR and R with strong stability. Conclusions: These results suggest that the new hybrid artifact detection and correction method enhances the viability of fNIRS as a functional neuroimaging modality.

and hemodynamic oscillations, which is generally known as neurovascular coupling. [1][2][3] As an inexpensive and portable optical neuroimaging technique, functional near-infrared spectroscopy (fNIRS) can noninvasively monitor brain function by simultaneously measuring the concentration of changes of cerebral oxyhemoglobin (Δ½HbO 2 ) and deoxyhemoglobin (Δ½Hb). [4][5][6][7] It has attracted increasing attention during recent years. A major advantage of fNIRS over neuroimaging techniques, such as functional magnetic resonance imaging and positron emission tomography, is that fNIRS is portable. Therefore, fNIRS technique is easily applied to long-term measurements, which are also much more likely to exhibit frequent movement and produce motion artifacts in comparison with the short-term experiments. For example, an increased prevalence of movement for long-term monitoring can induce artifacts in fNIRS data. Rapid head shaking introduces great amplitude and high-frequency fluctuation in fNIRS signals, whereas slow head rotation will cause slow sustained varying oscillation along with lasting baseline shift (BS). 8 Failure to correct for the artifacts may lead to biased or spurious conclusions. 7 Therefore, the motion artifacts are still problematic for fNIRS monitoring, particularly in longterm recordings with additional challenges to be addressed in comparison with the short-term experiments. And it is necessary to apply preprocessing approaches to eliminate the artifacts for fNIRS signal collection.
Some artifact correction approaches have been proposed for fNIRS measurements. However, each suffers from some drawbacks. For example, the simplest way is to discard the data segments polluted by the artifacts. 9 This approach is only suitable if the number of motion artifacts detected is low and the number of trials is high. However, it may not be ideal in cases where there are only a small number of trials. Adaptive filtering can be used for motion correction if the external signals are highly sensitive to motion artifacts but not to brain activity, such as an accelerometer. The drawback is that this approach cannot completely remove motion artifacts if the motion information in the external signals is not exact. [10][11][12] Robertson et al. compared different methods including adaptive filtering, wavelet filtering, ICA, 13 and linear regression, for removing motion artifacts from the measured signal, and claimed that the best motion artifact removal results were produced by ICA and regression. However, all these methods made use of the colocated channels to identify the signal that was related to motion. 14 Furthermore, the waveletbased (WB) methods could deal with motion spikes but exacerbate BS artifacts. 15,16 Correlationbased signal improvement can effectively remove large spikes caused by head motion and improve signal quality and spatial specificity. The method assumes that oxy-Hb and deoxy-Hb are perfectly negatively correlated and become positively correlated only when there is an artifact, although this assumption might not always be true during certain developmental stages or under abnormal brain physiology. 7 Metz et al. 17 presented an acceleration-based artifact reduction algorithm with data adaptive threshold specifically designed for long-term measurements. Accelerometer-based motion artifact removal (ABAMAR) can be easily adapted to various monitoring fNIRS applications for correcting baseline motion artifacts, with a drawback of the application of global threshold for the standard deviation (SD) of the amplitude or concentration signals. 8 Frank et al. 6 reported a temporal derivative distribution repair procedure to discard BS and spike artifacts. However, this technique mainly corrects the low-frequency signal, and it is still necessary to develop an approach to deal with the high-frequency component in fNIRS signals. 6 A spline interpolation method was proposed to identify and correct the BSs and smoothing methods to correct the spikes. 5 However, it might be difficult to define a global threshold of the maximum amplitude change for finding BSs, and many parameters need to be adjusted manually, such as signal-to-noise ratio (SNR) and the degree of spline function. Scholkmann et al. developed a method based on moving SD and spline interpolation to reduce movement artifacts. However, a limitation is that the variances of the motion artifacts have to be greater than that of the physiological signal to ensure the correct reduction of the artifacts. 18 Yucel et al. applied an approach in EEG applications to the applications of specialized NIRS optical fibers using improved optode-scalp coupling through miniaturized optical fiber tips. A potential limitation is that the position of the collodion-fixed fibers might introduce a bias. 19 Brigadoi et al. compared the performance of different motion correction techniques on real data during a cognitive task, and revealed wavelet filtering is the most effective approach. 20 Principle component analysis, spline interpolation, wavelet analysis, and Kalman filtering approaches were also compared using the accuracy of the recovered and simulated hemodynamic response function. The results showed spline interpolation produced the greatest improvement in mean-squared error, and wavelet analysis generated the greatest increase in contrast-to-noise ratio. 21 Thus each algorithm has its own advantage. However, how to combine the advantage of these approaches is still an issue. Regarding this, Hu et al. 22 demonstrated that the combined use of wavelet and moving average was preferable to each method on its own, whereas the combination of wavelet with targeted principal component analysis did not outperform wavelet alone. 23 However, the moving average technique cannot be considered a motion correction method such as wavelet since it served a high-pass filter function to remove slow drifts. Di Lorenzo et al. proved that the combination of spline and wavelet outperformed its individual use using infant semi-simulated data. However, motion artifacts in infant data greatly differ from adults' both in amplitude and frequency of occurrence. 24 The artifact correction recommendations derived from infant and children data might not be generally optimal for all kinds of fNIRS data. Thus every approach mentioned above has its own advantages and weaknesses.
WB methods can provide good localization property resulting in isolation of abrupt changes for the motion artifacts from hemodynamic changes. They are very effective in removing highfrequency spikes as well as low-frequency and low-amplitude artifacts but not in removing BSs. 16 On the other hand, spline interpolation can properly correct BSs by modeling the motions and subtracting from the original signal; however, it cannot deal with high-frequency spikes. 5 An optimum artifact removal method should be capable of detecting the artifacts objectively and correcting different types the artifacts efficiently. Here, we proposed a hybrid approach in this study for artifact detection and correction for fNIRS measurements, gathering the strengths of the existing approaches. The method first detects the artifacts through fNIRS detection strategy. The artifacts are divided into three categories: BS, slight artifacts, and severe artifacts. Finally, the artifacts are corrected based on the categories by three main steps: severe artifact correction by cubic spline interpolation, BS removal by spline interpolation, and slight artifacts reduction by dual-threshold WB method. In order to compare the proposed approach with other algorithms, we detected and corrected the artifacts in the fNIRS data acquired during whole night sleep monitoring using different methods. The metrics, SNR and Pearson's correlation coefficient (R), were applied for comparison.

Hemodynamic Signals Extraction
In this study, the hemodynamic signals were calculated from the optical signals acquired by fNIRS multiparameter collectors. Δ½HbO 2 and Δ½Hb were estimated using the modified Lambert-Beer Law 25 and the equations can be expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 2 8 5 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 2 3 0 where L represents the measuring depth of the fNIRS probes, and coefficients α 1 , α 2 , β 1 , and β 2 in the equations vary with the wavelength of near-infrared light.

Hybrid Motion Artifact Detection and Correction
In our study, the artifacts in fNIRS were divided into two categories: oscillation and BS. The oscillation should be further classified into two kinds: slight oscillation and severe oscillation based on its intensity. The slight and severe oscillations in the fNIRS signals correspond to the different levels of subjects' head movement severity, whereas the BS is generated simultaneously because of the change of head position. The artifact detection and correction approach put forward in this paper consists of artifact detection, severe oscillation correction, BS correction, slight oscillation correction, and high-pass filter, which can be described as follows (Fig. 1).
• Motion artifact detection. fNIRS signals were used to identify different kinds of motion artifacts including oscillations and BS.
• Severe motion artifact processing. Moderate and heavy oscillations, and BS removal were realized by spline interpolation.
• Slight motion artifact removal. Slight oscillation and the remnant motion artifacts can be dislodged with dual-threshold WB method and high-pass filtering.

fNIRS-based detection
The first step of the proposed method is the artifact detection using fNIRS-based detection strategy. fNIRS-based detection strategy uses fNIRS signal to detect oscillations and BS. Given a specific signal xðnÞ with sampling frequency Fs NIRS , the two-side moving SD tðnÞ of the measured signals can be calculated as 18 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 2 5 1 where n ¼ k þ 1; k þ 2; : : : ; N − k. N is the length of the time series xðnÞ, and k ¼ 3Fs NIRS . W is an odd value of the sliding window's length, and W ¼ 2k þ 1. A data segment for 240 s was selected from a continuous fNIRS signal that lasted 360 s to demonstrate the generation process of the ultimate dynamic threshold. By sorting the moving SD series tðnÞ in ascending order, pðnÞ is obtained. Then T 1 is defined as maximum value of the first 30% members, T 2 as maximum value of the first 50%, and T 3 as maximum value of the first 70%. Next, the data in tðnÞ are replaced with T 3 if exceeding T 3 , and with T 1 if lower than T 1 ; a modified moving SD series is acquired. The modified moving SD series is filtered by a high-pass filter with a relatively low cutoff frequency of 0.02 Hz. If the data value in the filtered series is lower than T 2 , it is set as T 2 . Then by adding 3ðT 3 − T 1 Þ, the ultimate dynamic threshold AðnÞ is obtained as shown in Fig. 2. The basic criterion for oscillation detection is defined as Eq. (4). t k and A k are the corresponding data point in moving SD series tðnÞ and dynamic threshold series AðnÞ, respectively. "None" means no oscillation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 5 1 9 As for suspected severe oscillation, further confirmation is performed. The oscillation period of fNIRS signal is divided into segments according to the slope changes. The principle is to ensure that the junction point between segments is located at the trough, and each segment contains only one peak. Another data segment for 35 s was selected to illustrate the confirmation process of the suspected severe oscillations in Fig. 3. Then the fNIRS signal range (x max − x min ) of a certain segment is compared with the SD of the suspected oscillation area, and x max and x min are the maximum and minimum values of the certain segment, respectively. If ðx max − x min Þ ≥ 6 × SD, the suspected severe oscillation can be confirmed as severe oscillation, otherwise the segments of the suspected severe oscillation should be detected as slight oscillation as shown in Fig. 3.
The BS is detected using fNIRS signal based on the oscillation detection results. If the time duration of two nonoscillation segments before and after the oscillation segment in the original fNIRS signal xðnÞ is longer than a certain length, which is 5 s selected in our study, the mean values during the two nonoscillation segments are calculated as ave before and ave after . x max and x min are the maximum and minimum values of the oscillation segment, respectively. If jave before − ave after j > ðx max − x min Þ∕2, a BS event could be marked for the oscillation area.  The original moving SD series t ðnÞ, the maximum value T 1 , the modified moving SD series, the filtered moving SD series, and the ultimate dynamic threshold AðnÞ.

Acceleration-based detection
In order to compare with the fNIRS-based detection, the signal from accelerometers can be used as an indicator of oscillation. First, the integrated acceleration signal accðnÞ from the M channels of the accelerometer is obtained E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 6 8 0 a k ðnÞ is the signal series of acceleration in different directions. In this study, M ¼ 3 and a k ðnÞ can represent a x ðnÞ, a y ðnÞ, and a z ðnÞ, respectively. The accðnÞ is measured with g (9.81m∕s 2 ). Next, the moving SD series dðnÞ of accðnÞ is sorted in ascending order. And the maximum value of the first 30% of ordered dðnÞ is defined as T. Then the value in dðnÞ is compared with T. If the specific data in dðnÞ exceed T, then it is set as T; if lower than 0.5T, it is set as 0.5T. Then the modified moving SD series d 0 ðnÞ is acquired. Finally, the dynamic threshold BðnÞ is obtained by doubling the d 0 ðnÞ. By comparing the original acceleration signals' moving SD series dðnÞ with BðnÞ, the severe and slight oscillation can be detected according to the same criterion as in fNIRS-based detection. As for BS detection, it is the same as the procedures in fNIRS-based detection.

Severe oscillation correction using cubic spline interpolation
The cubic spline interpolation implemented in MATLAB was used to correct severe oscillation in this paper. The cubic smoothing spline function is represented as f, and the interpolation process is calculated as 5,18 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; where xðnÞ is the signal needed to be processed, N is the length of xðnÞ, and D 2 f denotes the second derivative of the function f. The default value for the weight vector wðnÞ is all-one series. p is the smoothing parameter in the range of [0,1]. In our study, p was set to 0.99 as suggested in the previous studies. 5,18 The cubic spline interpolation was applied to the severe oscillation area to obtain s I ðnÞ. After subtracting s I ðnÞ from xðnÞ, x 1 ðnÞ was acquired.

Baseline shift removal using spline interpolation
As for BS removal, spline interpolation method was used with a different p ¼ 1∕ðFs NIRS Þ 3 . Fs NIRS is the sampling rate of the fNIRS signal. s S ðnÞ is the signal after applying spline interpolation to the BS area. After subtracting s s ðnÞ from x 1 ðnÞ and adding the beginning DC offset value within the BS area, the signal x 2 ðnÞ after severe artifact correction was obtained.

Slight oscillation correction using dual-threshold wavelet-based method
Slight oscillation and the remnant artifacts can be dislodged with dual-threshold WB method.
As isolated large wavelet coefficients in the discrete wavelet domain, slight artifacts can be easily identified and extracted by discrete wavelet decomposition. 16 Applying discrete wavelet transform (DWT) to x 2 ðnÞ, wavelet coefficients can be obtained E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 1 5 9 w jk ¼ w Ã jk þ n jk ; j¼ 1;2; : : : ; M: w Ã jk is the wavelet coefficient of the interested signal. n jk is the wavelet coefficient of the artifacts. A detailed procedure can be found in Ref. 16. The probability of is subordinated to normal distribution, and w jk can be described by a mixture of Gaussians functions. That is to say, w jk ∼ Nð0; σ 2 Þ and σ represents the SD of the mixed Gaussian distribution. All the w jk values are sorted in ascending order and numbered from 1 to L, and L represents the length of w jk . The data point whose number is nearest to 0.1587L can be found and marked asσ to function as the estimation of SD of all the w jk .
The artifacts usually result in outliers of wavelet coefficients through DWT. The artifact correction performance was compared using 16 commonly used mother wavelets. In this study, two key parameters w 1 and w 2 were proposed to determine where outliers are, which are called dual thresholds. w 1 is called detecting threshold to detect the artifacts existence. w 2 is named as processing threshold to determine whether the wavelets coefficients should be abandoned. The probability parameter α was used as a tuning parameter to determine the values of w 1 and w 2 . From Eq.
Φ is the normal cumulative distribution function. In this study, α was set from 0.005 to 0.1 in order to explore the influence of α on the effectiveness of the proposed method. w 2 ¼ u Ã jσj and w 1 ¼ u Ã 1.3jσj. Then the absolute value of w jk (jw jk j) is compared with w 1 according to its original order. If jw jk j is larger than w 1 , all the near wavelet coefficients whose absolute values are larger than w 2 are found from both sides and regarded as outliers. It is stopped when meeting the first wavelet coefficient whose absolute value is lower than w 2 . Then the values of all the outliers of wavelet coefficients are set to zero, and the target signal is obtained by wavelet reconstruction using processed wavelet coefficients (Fig. 4).
During the process of wavelet reconstruction, a cyclic averaging method is applied to avoid pseudo-Gibbs phenomena near singularities or abrupt changes and reduce the spike oscillation caused by the DWT as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 2 0 7 y 2 ðnÞ is the fNIRS signal after removal of severe oscillation, S is a shift operation, T represents the decomposition and reconstruction of discrete wavelets, and K is the number of cycles. In this paper, we selected K ¼ 16. The signals were filtered by a high-pass filter with the cut-off frequency of 0.003 Hz to remove the low-frequency oscillations. All the above algorithms used in our study were implemented in MATLAB and freely available upon request. Fig. 4 The processing of wavelet coefficients with dual-threshold WB method. The red circle areas indicate the outliers of wavelet coefficients exceeding the processing threshold should be processed. The blue triangles represent the coefficients between the two thresholds but it should not be changed because it has no neighbor whose absolute value is larger than detecting threshold.

Method Comparison
To validate the effects of the proposed approach, we compared the artifact removal performance of WB, ABAMAR, spline interpolation, median filtering, Spline Savitzky-Golay (Spline SG), spline-Rloess, and severe oscillation processed with cubic spline interpolation combined with wavelet filtering, and the proposed approach including artifact detection and correction. The ABAMAR, WB, spline interpolation, and cubic spline interpolation approaches were described as above. A one-dimensional median filter (window width = 600, step = 10) was applied to the raw temporal data to eliminate outliers such as those with sudden jumps or drops and remove impulsive noise and motion artifacts. The Spline-SG method is a multistep correction procedure. 5 A Sobel filter was first used to extract the gradient of the signal, from which time periods including motion artifacts are identified as those exceeding 1.5 times the interquartile interval. Then the BSs and slow spikes were identified and corrected those time periods using the spline approach. Finally, a Savitzky-Golay smoothing filter was performed to remove fast spikes. 5 As for the spline-Rloess technique, the identified BSs and slow motion artifacts were also corrected by the spline approach, and the spikes were removed using Rloess filters for the spline-Rloess method. The details can be found in Ref. 27.

fNIRS Dataset
In this study, a dataset was generated to evaluate the performance of the proposed approach. The study was conducted in accordance with the guidelines of declaration of Helsinki. All participants gave written informed consent. The fNIRS data were collected in the Sleep Disorders Diagnosis Center of Xijing Hospital, Military Medical University of The Air Force using a fNIRS measurement equipment developed by our laboratory to monitor the hemodynamic changes of the brain tissues at different depths in real time as shown in Fig. 5(a). 28 The details related to the equipment can be found in Ref. 28. The near-infrared probes of the equipment were symmetrically attached on the left and right forehead area of the subjects as shown in Fig. 5(b). The near-infrared signals were converted into the blood oxygen signals including surface deoxygenated hemoglobin (deoxy-Hb), surface oxygenated hemoglobin (oxy-Hb), proximal deoxy-Hb, proximal oxy-Hb, distal deoxy-Hb, distal oxy-Hb with units of μmol∕L, and local tissue blood oxygen saturation with a unit of percentage. An accelerometer was attached to each fNIRS probe to record the acceleration. The sampling frequency was 10 Hz. Forty healthy participants (16 male and 24 females, mean age 32 years) were involved in this study. The subjects were asked to sleep the whole night (approximately from 11:00 pm to 6:30 am) as usual until they were awoken by the researcher monitoring the measurements. Involuntary movements of the subjects during the whole night were not restricted.
According to the accelerometer's signals, a segment (totally lasted 1.5 h) of relatively clean distal signals (Δ½HbO 2 ) was selected as a criterion for the performance evaluation. Then 60 different segments (each lasted 1.5 h) of noise polluted Δ½HbO 2 signals were chosen from all the subjects by visual interpretation according to the accelerometer's signal recordings, in order to extract the noise signals. The accelerometer-based method 8 combined with WB method 16 and a 0.002-Hz high-frequency filter were applied to remove the high-frequency spike-like artifacts, low-frequency artifacts, and also the BSs in the signal. The processed signals were subtracted from the original signals to acquire the noise signals. Finally, 60 simulated segments of Δ½HbO 2 signals, each lasted 1.5 h, were generated by combining the selected clean signal with the 60 noise segments. A sample simulated time course without the artifact, with uncorrected artifacts, and corresponding corrected sample time course without the artifacts detection strategy, with the acceleration-based detection strategy, and with the fNIRS-based detection strategy were compared to validate the effects of artifact detection strategy. Finally, the comparison methods including WB, ABAMAR, spline interpolation, median filtering, Spline SG, spline-Rloess, and severe oscillation processed with cubic spline interpolation combined with wavelet filtering, and the proposed approach were applied to all of these data.

Evaluation Metrics
Two metrics, SNR and Pearson's correlation coefficient (R), were calculated to evaluate the effect of the artifact detection and correction method. The SNRs for the simulated noised signal and the processed signal after artifact correction are defined as SNR noised and SNR processed , respectively. Assuming that the clean signal is xðnÞ; n ¼ 1; : : : ; N, and XðnÞ is the noised signal, NðnÞ ¼ XðnÞ − xðnÞ, where NðnÞ is the noise signal. The SNR noised for the simulated noised signal XðnÞ is defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 4 8 2 SNR noised ¼ 10 log 10 P jxðnÞj 2 P jXðnÞ − xðnÞj 2 : As for the processed signals, yðnÞ is the processed signal after the artifact removal, and NðnÞ ¼ yðnÞ − xðnÞ. The SNR processed for the processed signal yðnÞ is defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 4 1 0 SNR processed ¼ 10 log 10 P jxðnÞj 2 P jyðnÞ − xðnÞj 2 : The R represents the similarity between the clean signal xðnÞ and the processed signal yðnÞ, which is calculated as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 3 4 0 3 Results

Effects of Probability Parameter and Mother Wavelet for Dual-Threshold Wavelet-Based Method
The effect of probability parameter α ranging from 0.005 to 0.1 was assessed for dual-threshold WB method on the performance of the proposed method as shown in Fig. 6. It indicates that the SNR and R value is almost the same, whereas α changes from 0.005 to 0.1. There is no significant difference on the SNR and R for α value from 0.005 to 0.1. The artifact correction performance was compared using 16 commonly used mother wavelets for dual-threshold WB method. The mother wavelet types from number 1 to 16 are shown in Table 1. Figure 7 demonstrates that only one type of mother wavelet (rbio3.3) performed relatively poor. There are significant differences on the SNR and R between rbio3.3 and other mother wavelets, whereas there is no significant difference for the other 15 types of mother wavelets. Figure 8 shows typical results of applying the dual-threshold WB method and wavelet filtering to a synthesized signal, respectively. The synthesized signal contains high-frequency    spike-like artifacts, BSs, and the amplitude changes with time. The dual-threshold WB method has the capability to remove those artifacts, and the wavelet filtering is poor for BSs.

Artifact Detection Strategies Comparison
A typical simulated wave and the corresponding waves obtained after severe oscillation and BS correction, and slight oscillation correction is shown in Fig. 9. A sample simulated time course without the artifacts, with uncorrected artifacts, and corresponding corrected sample time course without the artifact detection strategy, with the acceleration-based detection strategy, with the fNIRS-based detection strategy, and with the combination of acceleration-and fNIRS-based detection strategy are shown in Fig. 10, respectively. For the corrected sample time course without the artifact detection strategy, severe motion artifacts were not processed and only WB method was applied. Compared with the noised signal, the quality of the corrected signal was significantly improved through the artifact correction. Taking the clean signal as a reference, the corrected signal with the fNIRS-based detection strategy and with the combined strategy behaved best, nearly as good as clean data in contrast with the corrected signal without the artifact detection strategy and with the acceleration-based detection strategy. There is no significant difference between the corrected signals with the fNIRS-based strategy and the combined strategy (P > 0.05). Thus the fNIRS-based detection strategy is selected in our study. The simulated noised signal was compared with the corrected signal without detection strategy, with the acceleration-based detection strategy, with the fNIRS-based detection strategy, and also with the combined detection strategy, respectively. And the corrected signal without detection strategy was also compared with the signal with the acceleration-based detection strategy, with the fNIRS-based detection strategy, and with combined detection strategy, respectively.  The statistical analysis results across all the 60 fNIRS segments by separate paired t-test with Bonferroni correction are shown in Fig. 11. The mean and SD of the SNR and R for the corrected signal with the fNIRS-based detection strategy and with the combined detection strategy are almost the same. There are significant differences on the SNR and R between the corrected signal without detection strategy and with fNIRS-based detection strategy and also with the combined detection strategy (P < 0.05∕7).

Artifact Correction Methods Comparison Results
The artifact correction performance of the proposed method compared with different approaches, including WB, ABAMAR, spline interpolation, median filtering, Spline SG, spline-Rloess, and severe oscillation processed with cubic spline interpolation combined with wavelet filtering, are shown in Fig. 12. We could see the proposed method performed best both in the SNR and R. Statistical analysis was performed between the proposed method and other approaches including WB, ABAMAR, spline interpolation, median filtering, Spline SG, spline-Rloess, and severe oscillation processed with cubic spline interpolation combined with wavelet filtering, respectively. Separate paired t-test with Bonferroni correction for multiple comparison was applied to compare the proposed method with other approaches. The results proved there are significant differences on the SNR (WB, spline interpolation, and median filtering: P < 0.05∕7) and R (WB, spline interpolation, and median filtering: P < 0.05∕7) between the proposed method and other approaches. Fig. 11 The artifact correction performance comparison evaluated by SNR and R between the uncorrected signal, corrected signal without the artifact detection strategy, corrected signal with the acceleration-based detection strategy, and the corrected signal with fNIRS-based detection strategy. "o" represents the basic method; * means P < 0.05∕7. Fig. 12 The performance comparison of the proposed method with different approaches, including WB, ABAMAR, spline interpolation, median filtering, Spline SG, spline-Rloess, and severe oscillation processed with cubic spline interpolation combined with wavelet filtering, by (a) SNR and (b) R.
In single-threshold WB method, the value of the threshold will directly influence the number of outliers, resulting in the differences in the artifact correction effect. 15,16 However, a constant threshold value might not satisfy the requirements of all the signals. In this study, we analyzed the influence of the threshold in dual-threshold WB methods by comparing the results for different values of the threshold from 0.005 to 0.1. It can be found in Fig. 6 that the influence of the different on the performance of the proposed method is not significant both in SNR and R. In fact, the effect of varying from 0.005 to 0.2 on SNR and R was tested (the result is not shown in this paper). The results also show that proposed method is not sensitive to the change of. So the threshold changes will not directly influence the artifact correction effects for the dual-threshold WB method proposed in our study. The reason might be the interval between the detecting threshold and processing threshold makes the result of the artifactWe compared artifacts removal correction is not so sensitive to the threshold changes for the signals with different qualities. Thus the disadvantages of the single-threshold WB method can be overcome in the dual-threshold approach. It can be inferred that the dual-threshold approach can locate the artifacts exactly for the fNIRS signal with motion artifacts and avoid useful information be abandoned to some extent for the relatively clean signals.
Mother wavelet also performs an important role in DWT. We compared artifact removal performance of 16 different mother wavelets in our study as shown in Fig. 7. In order to keep the valid frequency band of interest (0.003 to 0.04 Hz) 29,30 and discard the low-frequency oscillation, which is the offset after BS removal, a high-pass filter with cutoff frequency at 0.003 Hz was used. We found that among 16 commonly used mother wavelets only one (rbio3.3) performed relatively poor. It revealed that the wavelet denoising effect of the proposed method in our study is not sensitive to the type of the mother wavelets.
We analyzed the performance of artifact removal with different artifact detection strategies. Taking the clean signal as a reference, it can be concluded that the corrected signal with the fNIRS-based detection strategy and with the combined strategy performed best, almost as good as clean data compared with the corrected signal without the artifact detection strategy and with the acceleration-based detection strategy as shown in Fig. 11. There is no significant difference between the corrected signals with the fNIRS-based strategy and the combined strategy. The acceleration-based detection strategy performed relatively poorer because not every fluctuation in acceleration signals corresponded to an oscillation or BS in the same position of fNIRS signals. So the acceleration-based detection strategy may result in a loss of some useful information. Similarly, not every fluctuation in fNIRS signals was caused by the motion. To make it simpler, we recommend detecting the artifacts based on fNIRS signals. Thus it could work even if there are no acceleration signals.
As we can see from Fig. 12, the proposed method performed better than WB, ABAMAR, spline interpolation, median filtering, spline, spline-SG, spline-Rloess, and severe oscillation processed with cubic spline interpolation combined with wavelet filtering method both in the SNR and R. Although Brigadoi et al. 20 reported that WB method was the most effective method to remove the low-frequency and low-amplitude artifacts compared with other methods, some researchers found that it might be inability to correct lasting shift when applied to fNIRS data, which contain various artifacts especially BSs. 15 We can see from Fig. 10 that WB method performance was worst among all the techniques. It might be caused by the fact that WB method cannot deal with large BSs. The application of dynamic threshold and segmentation in the artifact detection in our proposed method can locate the different types of the artifacts more accurately, then the artifact correction algorithm can be adopted to match the artifacts with different types. Furthermore, WB method using the median value of wavelet coefficients to estimate the SD is not reasonable especially when the signal's quality was extremely poor. WB method imposed a single Gaussian distribution on wavelet coefficients. As for the dual-threshold WB method in the proposed method, we improved the way to estimate the SD of the Gaussian distribution and introduced the concept of dual threshold. The distribution of wavelet coefficients can be described by a mixture of Gaussians. 31,32 Additionally, the effectiveness of WB method is sensitive to the threshold value. 15 In contrast, the dual threshold in the proposed method is appropriate to fNIRS signals in different qualities, and it can perform well with different thresholds as mentioned above. The ABAMAR performs well at the correction of BS, but it uses the mean value in place of the whole data segments with the artifacts, which ignores some valid information. The use of spline interpolation method can overcome the problem, which could correct the artifacts as well as reserve the useful information. The drawback of this approach occurs if the artifacts are difficult to detect. 20 The performance of spline-SG, spline-Rloess, and severe oscillation processed with cubic spline interpolation combined with wavelet filtering method is similar to the ABAMAR both in the SNR and R in our study. The Rloess and SG methods could remove high-frequency motion artifacts with relatively weak performance in correcting slow motion artifacts or BSs. 5 The proposed approach demonstrated effective performance on the removal of motion artifacts and BS. On the other hand, the noises and artifacts in fNIRS measurements are very challenging, and they might be from various sources with distinct statistical properties. For example, instrument noise, such as electronic (Johnson-Nyquist) noise and shot noise, may be induced by hardware and surrounding space with a Poisson distribution. 33 Low-frequency drift may be introduced by small instabilities in the laser diode light sources. Measurement noise is assumed to be Gaussian white noise. Finally, signal variation introduced by the human subject may include motion artifact and "global" oscillations from systemic sources (pulse, respiration, blood pressure Mayer waves, and other low-frequency oscillations), as well as the localized hemodynamic response to task-induced brain activation. Since photons reflected from brain tissue must travel through the skin, skull, and dura before reaching a detector on the scalp surface, there is also oscillation in the fNIRS signal arising from fluctuations in the superficial layers of the scalp. 14 Therefore, it is hard to say that the artifacts can be totally removed by pure signal processing.
In conclusion, considering the disadvantages of the previous approaches, the proposed approach provides a distinct detection of the artifacts through fNIRS detection strategy, and a comprehensive correction through three main steps: severe artifact correction by cubic spline interpolation, BS removal by spline interpolation, and slight oscillation reduction by dual-threshold WB method.

Disclosures
This study was approved by Ethical Committee of Xijing Hospital in Xian, China. The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.